! ----------------------------
!
! Main program.
!
! ----------------------------
program controllingTurbulence
  implicit none
	
  real*8 q
  real*8 k0 
  real*8 nu 
  real*8 dt
  integer*8 m

  parameter (q = 2)
  parameter (m = 22)
  parameter (dt = .12d-4)
  parameter (nu = 1.d-7)
  parameter (k0 = 2.d0**(-4))

  complex*16, dimension(1:m):: FORCE
  real*8,     dimension(1:m):: K, K2, A, B, C, GAMA
  
  call initForce(FORCE, m)
  call precalculate(K, K2, A, B, C, GAMA, nu, k0, q, m)
  call shellModels(GAMA, FORCE, A, B, C, dt, m)  
end program controllingTurbulence


! ----------------------------
!
! Computes the energy evolution with the given parameters and outputs this in 
! the file 'velocity.txt'.
!
! ----------------------------
subroutine shellModels(GAMA, FORCE, A, B, C, dt, m)
  implicit none 

  integer*8 i
  integer*8 j
  integer*8 n
  integer*8 m
  integer*8 iterBTS
  integer*8 iterOutput
  integer*8 numberLastBTS
  integer*8 numberTransientBTS
  real*8 dt
  real*8 dtol
  real*8 time
  real*8 outputTime
  real*8 bigTimeStep
  real*8,     dimension( 1:m  ):: A, B, C, GAMA
  complex*16, dimension( 1:m  ):: FORCE
  complex*16, dimension(-1:m+2):: U
  complex*16, dimension(-1:m+2):: F0, F1, F2, F3, F4, F5
  
  parameter (dtol = 1.d-8)
  parameter (outputTime = 0.1d0)
  parameter (numberLastBTS      = 4)
  parameter (numberTransientBTS = 20)
  parameter (bigTimeStep = 1000.D0)
  
  iterBTS = IDINT(bigTimeStep / dt / 6.d0 + dtol)
  iterOutput = IDINT(outputTime / dt / 6.d0 + dtol)
  print *, "Number of iterations in a Big Time Step = ",  iterBTS
  
  
  call initU(U, m)
  !call initAdamBashforth4(U, F0, F1, F2, F3, A, B, C, GAMA, FORCE, dt, m)
  call initAdamBashforth6(U, F0, F1, F2, F3, F4, F5, A, B, C, GAMA, FORCE, &
                          dt, m)

  do i = 1, numberTransientBTS
    do j = 1, iterBTS
      !call rungeKutta4(U, GAMA, FORCE, A, B, C, dt, m)
      !call rungeKutta6(U, GAMA, FORCE, A, B, C, dt, m)
      !call AdamBashforth4(U, F0, F1, F2, F3, A, B, C, GAMA, FORCE, dt, m)
      call AdamBashforth6(U, F0, F1, F2, F3, F4, F5, A, B, C, GAMA, FORCE, &
                          dt, m)
      print *, U(22)
    end do
    
    time = i * bigTimeStep
    print *, i, time
  end do

  open (10, file = 'very_long_velocity.txt')
  write (10, *) m
  write (10, *) numberLastBTS * iterBTS / iterOutput
  do i = 1, numberLastBTS * iterBTS / iterOutput
    do j = 1, iterOutput
      !call rungeKutta4(U, GAMA, FORCE, A, B, C, dt, m)
      !call rungeKutta6(U, GAMA, FORCE, A, B, C, dt, m)
      !call AdamBashforth4(U, F0, F1, F2, F3, A, B, C, GAMA, FORCE, dt, m)
      call AdamBashforth6(U, F0, F1, F2, F3, F4, F5, A, B, C, GAMA, FORCE, &
                          dt, m)
    end do
    
    write (10, *) time + i * outputTime

    do n = 1, m
      write (10, *) REAL(U(n))
      write (10, *) AIMAG(U(n))
    end do
  end do  
  close (10)
  
  return
end subroutine shellModels


! ----------------------------
!
! Initiates the vector U with proper values.
!
! ----------------------------
subroutine initU(U, m)
  implicit none

  integer*8 m
  complex*16, dimension(-1:m+2):: U

  U    = (0.d0,    0.d0)
  U(5) = (0.03d0, -0.065d0)
  U(6) = (0.06d0, -0.09d0)
  U(7) = (0.06d0, -0.09d0)
	
  return
end subroutine initU


! ----------------------------
!
! Initiates the vector FORCE with proper values.
!
! ----------------------------
subroutine initForce(FORCE, m)
  implicit none
  integer*8 m
  complex*16, dimension(1:m):: FORCE
    
  FORCE = (0.d0,0.d0)
  FORCE(4) = 5.d0 * (1.d0,1.d0) * 1.d-3
  
  return
end subroutine initForce


! ----------------------------
!
! Precalculates the vectors K, K2, A, B, C, GAMA
! using the arguments nu, k0, q, dt and m.
!
! ----------------------------
subroutine precalculate(K, K2, A, B, C, GAMA, nu, k0, q, m)
  implicit none

  integer*8 m
  integer*8 n
  real*8 q
  real*8 k0 
  real*8 nu
  real*8, dimension(1:m):: K, K2, A, B, C, GAMA

  do n = 1, m
    K(n)  = k0 * q**n
  end do

  K2   = K**2
  A    = K
  B    = -0.5d0 * K / q
  C    = -0.5d0 * K / q**2
  GAMA = nu * K2
  
  return
end subroutine precalculate
